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We perform fully general relativist ic (GR) simulations to address the inspiral and merger of 
binary white dwarf-neutron stars. The initial binary is in a circular orbit at the Roche critical 
separation. The goal is to determine the ultimate fate of such systems. We focus on binaries 
whose total mass exceeds the maximum mass (Mmax) a cold, degenerate EOS can support against 
gravitational collapse. The time and length scales span many orders of magnitude, making fully 
general relativist ic hydrodynamic (GRHD) simulations computationally prohibitive. For this reason, 
we model the WD as a "pseudo-white dwarf" (pWD) as in our binary WDNS head-on collisions 
study 1]. Our GRHD simulations of a pWDNS system with a O.OSM© WD and a IAMq NS show 
that the merger remnant is a spinning Thorne-Zytkow-like Object (TZIO) surrounded by a massive 
disk. The final total rest mass exceeds Mmax, but the remnant does not collapse promptly. To assess 
whether the object will ultimately collapse after cooling, we introduce radiative thermal cooling. 
We first apply our cooling algorithm to TZlOs formed in pWDNS head-on collisions, and show that 
these objects collapse and form black holes on the cooling time scale, as expected. However, when 
we cool the spinning TZIO formed in the merger of a circular-orbit pWDNS binary, the remnant 
does not collapse, demonstrating that differential rotational support is sufficient to prevent collapse. 
Given that the final total mass exceeds Mmax for our cold EOS, magnetic fields and/or viscosity 
may redistribute angular momentum, ultimately leading to delayed collapse to a BH. We infer that 
the merger of realistic massive WDNS binaries likely will lead to the formation of spinning TZlOs 
that undergo delayed collapse. 

PACS numbers: 04.25.D-,04.25.dk,04.40.Dg 



I. INTRODUCTION 

During inspiral and merger, compact binaries emit a 
large flux of gravitational waves (GWs), making them 
among the most promising sources for GWs detectable 
by ground-based laser interferometers like LIGO Q, 
VIRGO 0,[55, GEO [6^, TAMA f^,^ and AIGO as 
well a s by proposed space-based interferometers such as 
LISA [i3 and DECIGO pT|. Extracting physical infor- 
mation about these binaries from their GWs may shed 
light on determining their ultimate fate, but requires 
careful modeling of these systems in full general relativity 
(see [m for review and references therein). Most effort 
to date has focused on modeling black hole-black hole 
(BHBH) binaries (see [l3] and references therein), and 
neutron star-neutron star (NSNS) binaries (see [14] for 
a review), with some recent work on black hole-neutron 
star (BHNS) binaries in full GR 

As a follow-up to our investigation of binary WDNS 
head-on collisions [1], in this work we perform fully gen- 
eral relativistic simulations of circular- orbit WDNS bi- 
naries through inspiral and merger. Throughout we call 
this the "inspiral case" to distinguish it from the "head- 
on" collision case. WDNS binaries are promising sources 
of low-frequency GWs (for LISA and DECIGO) and, as 
we argued in [34] , possibly also high-frequency GWs (for 
LIGO, VIRGO, GEO, TAMA and AIGO), if the remnant 
ultimately collapses to a black hole. 

Like NSNS binaries, WDNS binaries are known to ex- 
ist. In [34] we compiled tables of 20 observed WDNS 



binaries and their measured orbital properties. The 
NS masses in these systems range between 1.26 Mq 
and 2.08 M©, and their distribution is centered around 
1.5 Mq. On the other hand, the WD masses in these 
systems range between 0.125 Mq and 1.3 M©, and their 
distribution is centered around 0.6 M©. Eighteen of 
these observed WDNS binaries have total mass greater 
than 1.65 M©, 8 of which have a WD component with 
mass greater than 0.8 M©, and 5 have total mass greater 
than 2.2 M©. This is interesting because the expected 
Tolman-Oppenheimer-Volkoff (TOV) limiting mass for a 
cold, degenerate gas must be larger than 1.97 Mq [35[ 
and may reach 2.2 Mq [36l-[4^, depending on the equa- 
tion of state (EOS). One of the main goals of this work 
is to determine whether a WDNS merger remnant will 
undergo prompt collapse to a black hole. 

Population synthesis calculations [45|, [Igj show that 
there are about 2.2 X 10^ WDNS binaries in our Galaxy, 
and that they have a merger rate between 10~^ yr~^ and 
1.4 X 10~^ yr~^. Furthermore, these studies find that af- 
ter a year of integration, LISA-like interferometers should 
be able to detect 1 — 40 WDNS pre-merger binaries. Re- 
cent work by Thompson et al. [43 suggests that the lower 
limit on the merger rate of binary WDNSs in the Milky 
Way, at 95% confidence, is 2.5 x 10~^ yr~^. Thompson et 
al. also suggest that the merger rate in the local universe 
is ~ 0.5 — 1 X lO^Gpc"^ yr"^. Therefore, ignoring some 
uncertainties, all recent population synthesis calculations 
suggest that LISA-like projects should be able to detect 
a few WDNS pre-mergers per year. 
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A. Previous WDNS work 

In we focused on possible evolutionary scenarios 
for circular WDNS binaries that have inspiraled suffi- 
ciently close that they reach the termination point for 
equilibrium configurations. This is the Roche limit for 
WDNSs, at which point the WD fills its Roche lobe and 
may experience one of at least two possible fates: i) sta- 
ble mass transfer (SMT) from the WD across the inner 
Lagrange point onto the NS, or ii) tidal disruption of the 
WD by the NS via unstable mass transfer (UMT). 

Note that once an UMT binary reaches the critical 
Roche separation, further inspiral and merger is governed 
by tidal effects and hydrodynamical interactions and not 
GW emission. 

We also studied key parameters that determine 
whether a system will undergo SMT or UMT and found 
that, for a given NS mass, there exists a critical mass ra- 
tio ^crit ^0.5 that separates the UMT and SMT regimes. 
If the mass ratio q = Mwd/^ns of a WDNS system is 
such that q > g'crit, the WD quickly overfills its Roche 
lobe, and the binary will ultimately undergo UMT. In the 
opposite case, q < ^crit, the system will undergo SMT. 
We showed that a quasistationary treatment is adequate 
to follow the evolution of an SMT binary during this secu- 
lar phase and calculated the gravitational waveforms. We 
also pointed out that WDNS observations suggest that 
there are known candidates residing in both regimes. 

In the case of tidal disruption (UMT), by contrast, 
the system will evolve on a hydrodynamical (orbital) 
time scale. In this scenario the NS may plunge into 
the WD and spiral into the center of the star, forming a 
quasiequilibrium confi gura tion that resembles a Thorne- 
Zytkow object (TZO) [4^; alternatively, the NS may be 
the receptacle of massive debris from the disrupted WD. 
WDNS mergers may also give rise to gamma ray bursts 

mil. 

The ultimate fate of the merged WDNS depends on 
(1) the initial mass of the cold progenitor stars, (2) the 
degree of mass and angular momentum loss during the 
WD disruption and binary merger phases, (3) the angu- 
lar momentum profile of the WDNS remnant, and (4) 
the extent to which disrupted debris is heated by shocks 
and/or nuclear reactions as it settles onto the NS and 
forms an extended, massive mantle. These are issues 
that require a hydrodynamic simulation to resolve. Note 
that Newtonian work on binaries with a WD component 
has been performed analytically in [s^^ l5ll - [55| and via 
Newtonian hydrodynamic simulations in [5o-m]. How- 
ever, ascertaining whether or not the neutron star ulti- 
mately undergoes catastrophic collapse (either prompt or 
delayed) to a black hole requires that such a simulation 
be performed in full general relativity. In fact, even the 
final fate of the NS in the alternative scenario in which 
there is a long epoch of SMT may also lead to catas- 
trophic collapse, if the neutron star mass is close to the 
neutron star maximum mass. This scenario too will re- 
quire a general relativistic hydrodynamic simulation to 



track. 

In [1] we employed the Illinois adaptive mesh refine- 
ment (AMR) relativistic hydrodynamics code [i^, EH to 
perform the first simulations of these systems in full GR. 
In particular, we studied the head-on collision from rest 
at large separation of a massive WD and a NS. We fo- 
cused on compact objects whose total mass exceeds the 
maximum mass supportable by a cold EOS in order to 
determine the outcome of such collisions. 

The vast range of time and length scales involved in 
the WDNS problem make fully general relativistic sim- 
ulations extremely challenging. In [1] we demonstrated 
that the length scales span four orders of magnitude, as 
measured in neutron star radii, and that the associated 
time scales span six orders of magnitude in M, the to- 
tal system mass. Current numerical relativity techniques 
and available computational resources make such simula- 
tions prohibitive. For this reason, we tackled this prob- 
lem using a different strategy. 

In particular, we constructed a six-parameter piecewise 
polytropic EOS which mimics realistic NS EOSs while, 
at the same time, scales down the size of the WD. We 
call these scaled-down WDs "pseudo-WDs (pWDs)". We 
chose all of the piecewise EOSs in such a way that the 
maximum NS mass is 1.8 Mq [63], and the maximum WD 
mass is 1.43 M©, i.e., the Chandrasekhar mass. Further- 
more, we made sure these EOSs preserve the qualitative 
shape of the central density-mass curves as well as the 
number of stable and unstable NS and WD branches (see 
Figs. 1 and 2 in [\\). Moreover, the scaling is performed 
so that all the length-scale and time-scale inequalities of 
the realistic problem are left unchanged. For a given set 
of EOS parameters, a realistic WD has a counterpart 
pWD which has the same mass but is smaller in size. As 
a result, for every realistic WDNS system, we can con- 
struct a pWDNS counterpart which involves the same 
(reahstic) NS and the pWD counterpart of the WD. 

Using pWDs we performed a sequence of head-on sim- 
ulations in which the EOS is changed so that the pWDs 
have the same mass (0.98 Mq) but decreasing com- 
pactions, while the compaction and mass of the NS in- 
volved remains practically unchanged. More precisely, 
while keeping the masses of the binary components and 
the NS radius fixed, the pWD compaction was modified 
so that the pWD:NS radius ratio varied between 5:1 and 
20:1. We then scaled the results of our simulations to 
predict the outcome in the realistic case: 500:1. 

In addition to studying the effects of the pWD com- 
paction, we also studied the effects of NS mass. We con- 
sidered NSs with masses 1.4 M©, 1.5 M©, and 1.6 Mq. 

All head-on collision simulations that we performed 
showed that after the collision, 14%- 18% of the initial to- 
tal rest mass escapes to infinity. In all cases, the remnant 
rest mass exceeded the maximum rest mass that our cold 
EOS can support (1.92 M©), and no case led to prompt 
collapse to a black hole. This outcome arises because 
the final configurations become hot, due to shock heat- 
ing. All our cases settle into a spherical quasiequilibrium 
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configuration consisting of a cold NS core surrounded by 
a hot mantle. Hence, all remnants are Thorne-Zytkow- 
like Objects (TZlOs). Scaling our results to realistic WD 
compactions, we predict that a realistic head-on collision 
will form a quasiequilibrium TZIO. 

Although the head-on collision simulations appear to 
lead to a consistent result (the formation of a TZIO), 
these results cannot be used to predict the final fate of 
WDNS systems in circular orbit. On the one hand, one 
might expect that the remnant in the inspiraling case 
will collapse to a black hole, because shock heating is 
not likely to be as intense as in the head-on case. On 
the other hand, the large amount of angular momen- 
tum in the inspiraling binary case may work to pre- 
vent prompt collapse. Therefore, to predict whether the 
merged WDNS remnant will collapse, promptly or fol- 
lowing cooling, we need to perform fully general rela- 
tivistic simulations of WDNS binaries through inspiral 
and merger. 



B. Goals and objectives 

The purpose of the current work is threefold: 
a) We simulate the late inspiral and merger of a WDNS 
system consisting of a IAMq NS and a O.OSM© WD ini- 
tially in circular orbit and at the Roche limit. As in our 
head-on colhsion studies, we employ the pWD approxi- 
mation to make the computations feasible [6^ . The pWD 
approximation is useful for predicting the ultimate fate 
of a realistic WDNS merger using scaling. In particular, 
the collision velocity {vc) and the pre-shocked WD sound 
speed Cs both scale as ~ (M/i^wo)^^^- This implies that 
the Mach number {A4 = Vc/cg) is invariant under scal- 
ing of RwB and so is the degree of shock heating. The 
thermal energy, as well as the rotational kinetic energy 
(T) and the gravitational potential energy {W) all scale 
as ~ M^/i?wD, when the binary merges. Thus, T/|W| 
is also invariant under scaling of Rwb- These consider- 
ations simply mean that with respect to gravity the rel- 
ative importance of thermal and rotational support in a 
WDNS merger remnant is approximately invariant, when 
the masses of the binary components are fixed and the 
only quantity that changes is the WD radius. As a con- 
sequence, the results obtained when adopting pWDNS 
systems can be scaled up to realistic WDNS systems. 
Note that our compaction study in [1] confirms the above 
scaling with the Mach number. 

b) We introduce a radiative cooling prescription and 
modify our adiabatic simulations by allowing for cooling 
to determine whether the merger remnant will collapse 
without thermal support, if it fails to collapse promptly. 
Otherwise, angular momentum provides sufficient sup- 
port to prevent collapse. 

c) We allow cooling to occur in the TZlOs formed in 
our WDNS head-on collision simulations [1] to confirm 
that these remnants collapse to a black hole when the 
excess thermal energy is radiated away. In other words. 



we demonstrate that it is thermal pressure alone that 
prevents these objects from undergoing prompt collapse, 
since angular momentum support is completely absent in 
head-on collisions. Delayed collapse occurs on a cooling 
time scale in all cases, providing a consistency check on 
our cooling implementation. 



Our pWDNS merger calculations show that the inspi- 
ral remnant is a spinning TZIO which is surrounded by a 
massive, extended, hot disk. In contrast to our head-on 
collisions, we do not find any outfiows in the inspiral- 
ing case. Therefore, the final total mass is greater than 
the maximum mass supportable by our cold EOS and 
many nuclear EOSs. However, the remnant does not col- 
lapse promptly to a BH. We find that the remnant is 
both thermally and centrifugally supported. To deter- 
mine whether centrifugal forces alone can support the 
remnant we incorporate cooling and find that the object 
does not collapse to a black hole. Therefore, the extra 
support provided by rotation is sufficient for preventing 
the collapse. 



Even though the TZIO does not collapse after cooling, 
we expect delayed collapse ultimately because the final 
total rest mass (~ 2.5 Mq) is larger than the maximum 
possible mass supportable by our cold EOS (and many 
nuclear EOSs), even allowing for uniform rotation. (The 
maximum gravitational mass of a uniformly rotating star 
with our adopted EOS is ~ 2.1 Mq). We expect that col- 
lapse to a BH will take place after viscosity or magnetic 
fields redistribute the angular momentum, as in the case 
of a hypermassive neutron star [65l-[67|. his conclusion 
will be true in the case of realistic WDNS mergers, un- 
less the true nuclear EOS supports a uniformly rotating 
star with a rest mass exceeding the remnant mass. Many 
viable EOSs do not support rest masses as large as 2.5Mq 
[37,], the remnant rest mass in our simulations. 

This paper is organized as follows. In Sec. [IT] the 
pWD approximation and the EOS adopted in our simu- 
lations are briefiy reviewed. Section IIIII outlines the ini- 
tial data generation technique. Sec. IIVI summarizes the 
methods used for evolving the gravitational and matter 
fields. Sec. [Vl introduces our radiative cooling formal- 
ism, which is then applied to the TZIO remnants from 
our pWDNS head-on collision simulations in Sec. IVII 
We present the results of our fully relativistic hydrody- 
namic simulations of the binary pWDNS late inspiral and 
merger in Sec. IVIIi and turn on cooling in Sec. IVII CI In 
Sec. I Villi we discuss possible effects of nuclear reactions 
in realistic WDNS mergers and give estimates of realis- 
tic cooling and angular momentum redistribution time 
scales. Sec. IIXI concludes with a summary of the main 
findings. Throughout this work, geometrized units are 
adopted, where G = c = 1, unless otherwise specified. 
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II. EQUATION OF STATE 



A. Gravitational Field Equations 



We employ the following 6-parameter piecewise poly- 
tropic cold EOS: 



P 



l/ni ^ 

1^1 Po ^ Po ^ Pi 

l^2po ^ Pl < Po ^ P2 

1/^3 

1^3 Po , P0> P2 



(1) 



where P is the pressure, po is the rest-mass density 
and /^i, /^2, /^3, ^1, ^2, ^3, Pi^ p2 are the parameters of the 
EOS. The parameters in Eq. ^ are 8 in number, but 
continuity requires that the following conditions be true 



l/n2 — l/ni 
1^1 l^2pi 



l/ns — l/n2 /o\ 
1^2 = I^3p2 • (2) 



As a result, the adopted EOS has 6 free parameters 
/^3, ^1,^2,^3,^1, and p2. 

Because of its multiple parameters, this EOS gives us 
the freedom to capture the same characteristic curves and 
turning points on a TOV mass-central density plot as for 
a realistic cold-degenerate EOS (see [68]), as shown in 
Fig. 1 in [1]. The EOS exhibits both stable {dM/dpo,c > 
0) and unstable {dM/dpo^c < 0) branches for both WDs 
and NSs, as in the reahstic case. 

Furthermore, this EOS allows us to adjust the size of 
a pWD of any given mass, thereby shifting the pWD 
branch to smaller radii (see Fig. 2 in Q), while keeping 
the NS branch approximately unchanged. For more de- 
tails about our EOS and pWDs we refer the interested 
reader to [l|. 

In this work the EOS parameters correspond to the 
10:1 EOS we considered in [1]: 1^3 = 4993, Fi = 
1.515, r2 = 2.969, Fa = 0.714, log(pi/pnuc) = -2.268, 
log(p2/Pnuc) = 0.208, where all values are in geometrized 
units and F^ = 1 + 1/n^, pnuc = 1.485 x 10~^km~^. These 
parameters are chosen such that the ratio of the isotropic 
radius of a TOV 0.98 M© pWD to that of a TOV 1.5 M© 
NS is 10:1. In addition, the EOS has been constructed so 
that the maximum gravitational mass of a NS is 1.8 M©, 
i.e., the same as that for the AP2 version of the Akmal- 
Pandharipande-Ravenhall (APR) EOS [3i,[69j, and the 
maximum gravitational mass of a pWD is 1.43 M©, i.e., 
the maximum mass of a TOV WD obeying the Chan- 
drasekhar EOS for mean molecular weight /ig = 2. 



III. INITIAL DATA 

This section introduces the formalism adopted for gen- 
erating valid general relativistic initial data for binary 
pWDNS systems in circular orbit. 



The spacetime metric in the standard 3+1 form [70| is 
written as 

ds^ = -a^dt^ + -fij{dx' + p'dt){dx^ + p^dt), (3) 

where a is the lapse function, the shift vector and 7^^ 
the three-metric on spacelike hypersurfaces of constant 
time t. Throughout the paper Latin indices run from 1 
to 3, and Greek indices run from to 3. 

The three-metric 7^^ is then conformally decomposed 

as 



7,, =^%-, 



(4) 



where ^ is the conformal factor and fij the conformal 
metric. We adopt the standard approximation of a con- 
formally flat spacetime, so that fij = 5ij in Cartesian 
coordinates. 

We split the extrinsic curvature {K'^^) into trace (K) 
and tracefree parts (A^^) 



o 

take the initial slice to be maximal 
K = 0, 



(5) 



(6) 



and introduce a "conformal", traceless extrinsic curva- 
ture as 



(7) 



Using Eqs. ©-([T]) and assuming the existence of 
an approximate helical Killing vector, the Hamiltonian 
and momentum constraint equations assume the form of 
the conformal-thin-sandwich (GTS) equations [l2|. The 
Hamiltonian constraint becomes 

V^^ = --^-^AijA'^ - 27r^V. (8) 
8 

where is the flat Laplacian operator associated with 
fij. Here the source term p is defined as 



p = n^n^T^^ 



(9) 



where is the normal vector to a t = constant slice, 
and is the stress-energy tensor of the matter. 
The momentum constraint yields 

V2/3^ + = 2A^\/j{a^-^) + 167ra^'^/, (10) 



where the source term f is given by 



(11) 
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TABLE I. Outer boundary conditions imposed on the CTS 
variables when generating WDNS initial data. 



Variable 


Fall-off condition 


^ - 1 


1/r 


a — 1 


1/r 


QX 




QV 


^ x/r^ 




^ xyzjr^ 


B 


^ xy/r^ 



Taking the trace of the evolution equation for Kij (see 
Eq. (2.106) in [12]), imposing the maximal slicing condi- 
tion Eq. (j6j), and combining the result with Eq. (j8j), we 
obtain an equation for the lapse [f2| 



V\a^) = -^-""AijA'^ + 27r^'^(p + 25) . (12) 



Here the source term S is defined as 
In all equations above 



(13) 



= ^ h'p' + ^'p' - Ir^kP" ] , (14) 



and Aij = fimfjn^^^- 

Instead of solving Eq. ([TOj) for the shift vector directly, 
it is convenient to decompose as a sum of a vector and 
a gradient (cf. [7l|) 
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(15) 



Eq. ([TQ|) can then be replaced by the two equivalent equa- 
tions 



and 



(16) 



(17) 



Equations (g]), ([Hj), ^ and ([I7j) form a system of 6 
coupled, nonlinear elliptic equations for the 6 unknowns 
^, a^, and 5, which must be solved iteratively. These 
equations are elliptic and hence require outer boundary 
conditions to be specified. We impose the same fall- 
off boundary conditions as in [72^, except that here we 
choose the binary components to be initially lined up on 
the X-axis and the binary rotation axis parallel to the 
z-axis. Table H] lists the full set of outer boundary condi- 
tions imposed in our initial data. 



B. Matter fields 

As we argued in [34] the WD in a WDNS binary with 
close separation likely will be tidally locked. For this 
reason we focus on corotating WDNS systems only. 

We assume that the matter is described by a perfect 
fluid stress-energy tensor: 



= (po + Pi + P)u''u^ + Pg 



a(3 



(18) 



where g' 



a(3 



the inverse of the four-metric and 



pQ^ pi^P^u'^ are the rest-mass density, internal energy 
density, pressure, and four- velocity of the fluid respec- 
tively. For all initial configurations, the pressure is given 
by the cold EOS as specified in Eq. ([1]). The internal 
energy density can be derived by integrating 



Po 



(19) 



and for Eq. ([T]) the integration yields 

Po ^ Pi 
f C2, pi < Po ^ P2 
C3, Po > P2 



Pi_ 

Po 



i/m 



l/n2 
n2^2Po 



(20) 



l/ns 

nKpQ 



where 

/ \ 1/^1 

C2 = [ni - n2)f^iPi , 



C3 = C2 + (n2 - ns)tZ2P2 



l/n2 



(21) 

In Cartesian coordinates we choose the orbital plane 
of the binary to be the z = plane, so that the fluid 
four- velocity takes the form ^| 



u\l, -ny,Q{x - Xrot),0), 



(22) 



where Q is the constant orbital angular velocity and Xrn t 
is the X coordinate of the axis of rotation. Following [72|, 
we introduce a vector 

r = (0,-?/,(x-Xrot)), (23) 

and rewrite the four- velocity as 

ix" = u\an'' + + f^"")' (24) 

The source term p in Eq. ([9|) can then be written 

Po + Pi + ^ 



P. 



(25) 



where v is the magnitude of the three- velocity of the fluid. 
Using u^Ua = —1, it can be shown that v'^ is given by 



^4 



{Qy - p-f + {^l{x - xrot) + + (r )' 



(26) 
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The momentum source in Eq. (pTj) becomes 

, _ (po + ^, + p) {ne + 

^ " a l-t;2 ' ^^^^ 

and S in Eq. (fT^ is given by 

.5 = (/9o + Pi+P)-^ + 3P. (28) 
1 — 

C. Computational methods 

We solve the nonlinear elliptic equations (j8j), ([T2j), 
(p!6|) and ([T7j) using a fixed-mesh-refinement (FMR) fi- 
nite difference code we developed, which is based on the 
Portable, Extensible Toolkit for Scientific Computation 
(PETSc) library jTMEj. A fuh description of our code 
may be found in Here we summarize the basic fea- 
tures. 

The grid structure used in our FMR elliptic code is 
a multi-level set of properly nested, uniform grids. We 
use standard cell-centered, second-order accurate finite 
difference stencils for the Laplacian operator and the 
derivatives of the variables, using first-order interpola- 
tion across the refinement level boundaries when neces- 
sary. We calculate the solution across the entire grid, 
and only on leaf cells (i.e. cells within which there exist 
no higher resolution cells). In p]] we performed a series 
of tests involving single NSs, and we demonstrated that 
the code converges to the expected solutions at second 
order. 

Given the matter distribution, Q and Xrot we solve the 
CTS equations iteratively, addressing the non-linearity of 
Eq. (jSj) by performing Newton- Raphson iterations, until 
the residuals of all six equations become smaller than 
some set tolerance (usually set to 10~^^). 

We obtain the WD rest-mass density distribution, Q 
and Xrot at the Roche limit for equilibrium, corotating 
binary WDNSs in circular orbit obeying our cold EOS us- 
in g th e unigrid Newtonian code we developed and tested 
in [3^. At the Roche limit, the binary separation is large 
enough so that the tidal effects on the NS are negligible, 
and hence the NS will be spherical to a high degree and 
point-like from the point of view of the WD. Thus, in the 
Newtonian code we model the NS as a point mass and 
we self-consistently solve for the WD rest-mass density 
distribution via the integrated Euler equation. We use 
the Newtonian equations for this step, because it is com- 
putationally simple and fast. Also, the large separation 
at the Roche limit ensures that the WD and NS inter- 
action lies in the Newtonian regime, so that our initial 
configuration is nearly in equilibrium. 

After the WD rest-mass density distribution has been 
calculated, the point-mass NS is replaced by a TOV NS 
with gravitational mass equal to that of the point-mass 
NS, centered at the position of the point mass. For sim- 
plicity, we model the NS as corotational because there 
is no essential difference between an irrotational and a 



corotational NS at such large separations. The spin of 
a corotating NS is very small. To understand this, con- 
sider the ratio of the angular velocity of the corotating 
NS (l^cor) to that at the mass-shedding (l^ms) hmit: 

^cor / ^total f RnS \ -,^fRNs\^^^ /^^x 

where Ar is the Roche limit separation. For the typi- 
cal system we consider i?Ns/^R ^ ^Ns/3i^wD- For re- 
alistic massive WDs Rwd/Rns ^ 500, and for pWDs 
RpWb/Rns ^ 10- Thus, l^cor/^ms ^ 10~^ for realistic 
WDNSs and l^cor/^ms ^ 10"^ for pWDNSs. Therefore, 
the corotation spin the NS acquires is very small and has 
no physical significance. 

Having prescribed the NS and pWD rest-mass density, 
using second-order polynomial interpolation, we interpo- 
late the NS and pWD matter distribution on the grid of 
our FMR elliptic code and solve the CTS equations. 



IV. EVOLUTION OF WDNS SYSTEMS 
A. Basic Equations 

The formulation and numerical scheme for our simula- 
tions are the same as those reported in [l|, [H, [62|, |76| , to 
which the reader may refer for details. Here we introduce 
our notation and summarize our method. 

We use the 3+1 formulation of general relativity, in 
which the metric is decomposed as in Eq. ([3|). In this 
formalism, the fundamental dynamical variables for the 
metric evolution are the spatial three-metric and ex- 
trinsic curvature Kij. The Baum^arte- Shapiro- Shibata- 
Nakamura (BSSN) formalism [l2|,[73,|7^ is adopted. The 
BSSN evolution variables are the conformal exponent 
(/) = ln(7)/12, the conformal 3-metric 7^^ = e~^^^ij^ 
three auxiliary functions F* = — 7*-^,j, the trace of the 
extrinsic curvature and the trace- free part of the con- 
formal extrinsic curvature Aij = e~'^^{Kij — ^yijK/S). 
Here 7 = det (jij). The full spacetime metric g^i, is re- 
lated to the three- metric 7^^^ by j^jy = ^^r/ + n^njy, where 
the future-directed, timelike unit vector normal to the 
time slice can be written in terms of the lapse a and shift 
as = a~^(l, — The evolution equations of these 
BSSN variables are given by Eqs. (9)-(13) in @. 

We adopt standard puncture gauge conditions: an ad- 
vective "1+log" slicing condition for the lapse and a 
"Gamma- freezing" condition for the shift [t^. Thus, we 
have 

doa = -2aK , (30) 
dof3' = (3/4)5^ , (31) 
doB' = dor - rjB' , (32) 

where do = dt — f^^dj. We set the rj parameter to 
0.01 km~^ for all simulations presented in this work. 
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The fundamental matter variables are the rest-mass 
density po^ specific internal energy e, pressure P, and 
four- velocity u^. We write the stress-energy tensor as 

T^u = Pohu^Uj, + Pg^j, , (33) 

where = 1 + e+P/po is the specific enthalpy and e is the 
specific internal energy. In our numerical implementation 
of the hydrodynamics equations, we evolve the following 
"conservative" variables: 

P* -^/ipon^u^ , (34) 

~S, ^ -VlT,,n^Y^ , (35) 
f = ^T^,n^n^ - . (36) 

The evolution equations for these variables are given by 
Eqs. (27)-(29) in [62]. 

The EOS we adopt for the evolution has both a thermal 
and cold contribution, and can therefore be written 

P = Pth + Pcoid, (37) 

where Pcoid is given by Eq. ([T]) and the thermal pressure 
is given by 

Pth = (rth-l)/>o(e-ecoid), (38) 

where 

Ccoid = - j Pco\dd{l/ Pq). (39) 

We set Fth = 1-66 {c:^ 5/3) in all our simulations. That 
is, we set Fth to the Fi exponent of the 10:1 EOS, appro- 
priate either for nonrelativistic cold, degenerate electrons 
or (shock) heated, ideal nondegenerate baryons. Equa- 
tion (|37|) reduces to our piecewise polytropic law Eq. ([1]) 
for the initial (cold) NS and pWD matter. 



that employs the piecewise parabolic (PPM) reconstruc- 
tion scheme [82j, coupled to the Harten, Lax, and van 
Leer (HLL) approximate Riemman solver [83]. The 
adopted hydrodynamic scheme is second-order accurate. 
To stabilize our hydrodynamic scheme in regions where 
there is no matter, a tenuous atmosphere is maintained 
on our grid, with a density floor patm set to 10~^^ times 
the initial maximum density on our grid. The average 
density of the pWD is lO^^gr/cm^, and at least six or- 
ders of magnitude larger than that of the artificial atmo- 
sphere. Thus, the atmosphere poses no problem in evolv- 
ing the pWD. The initial atmospheric pressure Patm is set 
by using the cold EOS ([1]). Throughout the evolution, we 
impose limits on the pressure to prevent spurious heating 
and negative values of the internal energy e. Specifically, 
we require Pmin < P < Pmax, where P^ax = lOPcoid and 

Pmin 

= 0.8Pcoid, where Pcoid is the pressure calculated 
using the cold EOS ([1]). Whenever P exceeds Pmax or 
drops below Pminj we reset P to Pmax or Pmin, respec- 
tively. Following |2j we impose the upper pressure limits 
only in regions where the rest-mass density remains very 
low (po < lOOpatm), but wc imposc the lower limit every- 
where on our grid. We impose the pressure floor every- 
where, because numerical error sometimes leads e — ecoid 
slightly below zero, resulting in negative thermal pres- 
sure. We have found experimentally that if this situation 
arises, it can be avoided in the subsequent timesteps by 
imposing the pressure floor. 

At each timestep, the "primitive variables" po, and 
v'^ must be recovered from the "conservative" variables 
p*, f, and Si. We perform the inversion numeri cally as 
specified in [62^ . We use the same technique as in [sj [85| 
to ensure that the values of Si and f yield physically valid 
primitive variables. 



C. Diagnostics 



B. Evolution of the metric and hydrodynamics 

We evolve the BSSN equations using fourth-order ac- 
curate, cell-centered finite-differencing stencils, except on 
shift advection terms, where fourth-order accurate up- 
wind stencils are applied. We apply Sommerfeld out- 
going wave boundary conditions on all BSSN fields, as 
in [23|. Our code is embedded in the Cactus paral- 
lelization framework [80|, and our fourth-order Runge- 
Kutta timestepping is managed by the MoL (Method of 
Lines) thorn, with the Courant-Friedrichs-Lewy num- 
ber set to 0.45 in all pWDNS simulations. We use the 
Carpet [8l| infrastructure to implement the moving-box 
adaptive mesh refinement. In all AMR simulations pre- 
sented here, we use second-order temporal prolongation, 
coupled with fifth-order spatial prolongation, and impose 
equatorial symmetry to reduce the computational cost. 

We write the general relativistic hydrodynamics equa- 
tions in conservative form. They are evolved via a high- 
resolution shock-capturing (HRSC) technique [62|, [76| 



During the evolution, we monitor the normalized 
Hamiltonian and momentum constraints as defined in 
Eqs. (40)-(43) of [23]. We also monitor the ADM mass 
and angular momentum of the system. The equations 
used to calculate the ADM mass and angular momen- 
tum with minimal numerical noise are as follows p^] 

-^X^A umdnl + ^^X^Sn] • (41) 

Here = e^, p = n^rij.T^'', Si = -n^7i^T^^, R is the 
Ricci scalar associated with 7^^, and Tijk are Christoffel 
symbols associated with jij. 
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TABLE 11. Summary of initial configurations. Mns (Mwd) is the ADM mass of an isolated NS (pWD)^"^ i?NS (i^wo) the 
isotropic radius of an isolated NS (pWD), Cns the compaction of an isolated NS, where the compaction is the ratio of the 
ADM mass of the isolated star to its areal radius. All pWDs considered here have Cwd = 0.01. Madm is the ADM mass of the 
system and A the initial binary separation in isotropic coordinates. Cases Al and A3 are the same as the head-on collision cases 
we studied in 1]. The initial coordinate separation for these cases was set to 587 km. Case A corresponds to our simulation 
of a binary pWDNS in circular orbit starting at the Roche limit. For this case QMadm 6.95 X 10"^. All cases have been 
produced with the 10:1 EOS of [1]. 



Case 


Mns/M© 


Mwb/Mq 






Rwb/Rns 


i?WD/MADM 


Madm /Mo 


A/i?WD 


Al 


1.4 


0.98 


0.11 


0. 


8.88 


41.18 


2.41 


4.00 


A3 


1.6 


0.98 


0.15 


0. 


11.15 


37.46 


2.65 


4.00 


A 


1.4 


0.98 


0.11 


2.88 


8.88 


40.05 


2.48 


3.14 



^"-^ Here we list the ADM masses, isotropic radii and compactions of the isolated NS stars, whose rest-mass density profiles 
were used to generate initial data. The same holds for the pWDs in cases Al and A3. In case A the pWD rest-mass density 
profile, the Roche limit separation and Q were generated by a Newtonian binary WDNS code and then used in our CTS solver. 



When hydrodynamic matter is evolved on a fixed uni- 
form grid, our hydrodynamic scheme guarantees that the 
rest mass Mq is conserved to machine roundoff error. 
This strict conservation is no longer maintained in an 
AMR grid, where spatial and temporal prolongation is 
performed at the refinement boundaries. Hence, we also 
monitor the rest mass 

Mo = y p^d^x (42) 

during the evolution. Rest-mass conservation is also vio- 
lated whenever po is reset to the atmosphere value. This 
usually happens only in the very low-density atmosphere. 
The low-density regions do not affect rest-mass conser- 
vation significantly. 

In all simulations we present in this work the nor- 
malized Hamiltonian constraint violations remain smaller 
than 0.9% and the normalized momentum constraint vi- 
olations smaller than 2.1%. Rest mass is conserved to 
within 4% and angular momentum to within 10%. 

Shocks occur when the stars collide. We measure 
the entropy generated by shocks via the quantity K = 
P/Pco\d ^ 1, where Pcoid is the pressure associated with 
the cold EOS that characterizes the initial matter (see 
Eq. (HD). 



V. RADIATIVE COOLING 

Our binary WDNS head-on collision studies in [1] 
demonstrate that the hot, quasiequilibrium TZIO rem- 
nants do not collapse promptly to a black hole, even 
though the final total mass is larger than the maximum 
mass supportable by the cold EOS. This outcome might 
also arise in the case of WDNS mergers in circular or- 
bit, and can be due to additional support provided by 
thermal pressure and/or rapid rotation. In order to de- 
termine whether thermal support is dominant, we add 
radiative cooling to the GR hydrodynamic equations. We 
now describe our formalism for implementing this. 



The dynamics of radiation is governed by (scl-lssj 

VaR""^ = -G^. (43) 
where R^^ is the radiation stress-energy tensor given by 

i?"^ = j dvdVtl^N'^N^, (44) 

and is the radiation four-force density given by 

= dudn{xuIu-ju)N''. (45) 

In the equations above dQ is the solid angle, u and 
Iiy = Iu{x^ ^ N\u) are the radiation frequency and spe- 
cific intensity of radiation at x'^ moving in direction 
j\^a _ jfiu^ respectively. All quantities are measured 
in the local Lorentz frame of a fiducial observer with four- 
velocity ii^.^, i.e., 

hv = -Paujid^ (46) 

where is the photon four-momentum and h denotes 
Planck's constant. The energy- momentum conservation 
equation then becomes 

V«(T"^ + i^''^) = (47) 

or after using Eq. (|43|) 

Vc,T"^ = G^. (48) 

After projecting this equation using the fluid four- 
velocity ii", we obtain the modified energy equation: 

U^'WaS = -{e + P)\/aU'' - u'^Ga, (49) 

where the perfect fluid stress-energy tensor has been writ- 
ten as 

T^^ = (e + P)^x^ii^ + Pg^'^ (50) 

and £ is the total energy density. Using the continuity 
equation 

Vc.(po^") = 0, (51) 
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TABLE III. Grid configurations used in our simulations. Here M is the sum of tfie ADM masses of the isolated stars, Ax 
is the grid spacing in the innermost refinement box surrounding the NS, A^ns denotes the number of grid points covering the 
diameter of the NS initially, and A^vD denotes the number of grid points covering the (smallest) diameter of the pWD initially. 
The outer boundary distance to the center of mass is approximately 1020M in cases Al and A3, and 540M in case A. 



Case M/Mq 


Grid Hierarchy (in units of M)*^""^ Ax 


A^NS AwD 


Al 2.38 


(534.33, 267.16, 133.58, 66.79, 35.78, 19.08, 10.44, 7.156) M/6.71 


63 


35 


A3 2.58 


(467.27, 233.64,116.82, 58.41, 29.20 , 15.58, 8.518, 5.841) M/8.22 


56 


38 


A 2.38 


(270.87, 135.44, 67.72, 36.28[N/A], 19.35[N/A], 9.674[N/A], 5.744[N/A]) M/13.2 


124 


73 



^"'^ There are two sets of nested refinement boxes: one centered on the NS and one on the pWD. This column specifies the 
half-length of the sides of the refinement boxes centered on both the NS and pWD. When there is no corresponding pWD 
refinement box (as the pWD is much larger than the NS), we write [N/A] for that box. 



Eq. (j49j) becomes 



(52) 



The total energy density is related to the specific thermal 
energy via the following equation 



Po(l + eth + ecoid). 



(53) 



Using Eqs. (|52)) and ([53|) we find that the specific thermal 
energy evolves as 



1 



Po 



po 



(54) 



where we have used dccou/dpo = Pcou/ Po, and Pth = 

P-Pcold- 

In the comoving reference frame w"Va = d/dr, where 
r is the proper time. Thus, in the comoving frame Eq. 
becomes 



A. 

d^ 



Pth d 



Po 



-po 



Po 



(55) 



In order to achieve cooling, the radiation term in (|55]) 
must be specified so that thermal energy can be removed. 
For this reason we choose the following cooling law that 
gives rise to exponential cooling 



u^Ga = ethPo/'Tc, 



(56) 



where Tc > is the cooling time scale. Substituting 
Eq. ([56|) in Eq. (|55]) we obtain 



d 
dr 



(Fth - 1) dpQ 
dr 



po 



(57) 



where we used Eq. (|38]) to substitute for the thermal 
pressure. 

The first term in brackets on the RHS of Eq. (|57|) re- 
sults from adiabatic compression or expansion. The sec- 
ond term results from cooling and radiates away ther- 
mal energy exponentially. Thus, if initially we have a 



quasiequilibrium spherical object which is thermally sup- 
ported, and we cool it quasistatically, i.e., choose a cool- 
ing time scale much longer than the free-fall time scale 
of the star, it will radiate thermal energy away and con- 
tract. While the contraction generates extra heat, radia- 
tion tends to remove any residual thermal energy. Thus, 
after cooling, TZlOs are expected either to collapse to 
a BH or be supported by the residual cold pressure and 
centrifugal force. 

In the optically thin regime, assuming the fluid radi- 
ates isotropically in its rest frame, and using u%(j = 
the source term G^ is generally expressed as [83, L 



G" = M"(r - A) + y dHxl + xt)F^, (58) 

where F and A are the heating and cooling terms, respec- 
tively, given by 



r = j dud^xlL. 

A = / dvdVLjy^ 



(59) 



and where Xv->Xt ^i"^ the absorption and scattering co- 
efficients, respectively, and jjj is the emissivity. Finally, 
is the total radiation flux four- vector 



(60) 



where the projection operator P"^ is defined as 

P";3 = 5";3+m"u^. (61) 

If we assume that there is no absorption and no scat- 
tering, Eq. (|58|) becomes 



As a result, 



u^G"" = A. 



(62) 



(63) 
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FIG. 1. Snapshots of rest-mass density profiles at selected times for head-on case Al after cooling is turned on. The contours 
represent the rest-mass density in the orbital plane, plotted according to po = po,maxlO~°'^*''~°'^® (j = 0, 1, . . . , 7). The color 
sequence dark red, red, orange, yellow, green, light green, blue and light blue implies a sequence from higher to lower values. 
This roughly corresponds to darker grey-scaling for higher values. Here po.max = 4.645/5nuc, where pnuc = 2 x 10"''* g/cm^. 
These snapshots clearly demonstrate that the entire TZIO remnant collapses once cooling is turned on. Here M = 2.38 Mq = 
3.52 km = 1.17 x 10"*^ s is the sum of the ADM masses of the isolated stars. 



A straightforward comparison of Eqs. ([56|) and ([63]) shows 
that the integrated emissivity in our coohng model is 
given by 



A — eth. 

Tc 



(64) 



If we project Eq. (|48|) using the timelike unit vector 
normal to spacelike hypersurfaces and the projection 
operator h^f^ = 5^ f3^n^n(3^ we find that the 3+1 GRHD 
equations become 

dtS, + dj{a^TU) = \(^VlT''^9a^.^ - o^^/lu^K, (65) 



and 



dtf + di{a^^T^' - p,v') = s - a^^u^A, (66) 



where we have used Eq. (|62]) and A is given by Eq. (j64|) . 
Thus, cooling enters as a source term in the GRHD equa- 



tions, which is precisely how cooling is implemented in 
our HRSC code. 

Note that the optically-thin approximation employed 
here is valid for neutrino cooling in the WDNS merger 
scenario (head-on or otherwise). According to our analy- 
sis in [1], the temperatures and densities of the hot mantle 
of a TZIO are such that thermal neutrino emission likely 
will be the dominant source of cooling. The diffuse TZIO 
mantle composed of the WD debris is optically thin to 
neutrinos, justifying the above approximation. 

Finally, note that the self-gravity of the radiation is ne- 
glected, i.e., we assume that radiation does not affect the 
spacetime structure, so only the perfect fluid stress en- 
ergy tensor contributes to the BSSN source terms. This 
is a good approximation as long as the radiation energy 
density is subdominant (n^n/^i?"^ <C nan^T^^). This 
is indeed the case in a WDNS scenario because the rest 
mass dominates the mass-energy as can be inferred by 
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FIG. 2. Left panel: Evolution of maximum value of rest-mass density with cooling (solid curve) and without cooling 
(dotted curve) for the case shown in Fig. (TJ Here po,max is the maximum value of the rest-mass density, y9o,max(tc) = 5.879 x 
10"^ km"2 = 7.919 x 10^^ g/cm^ is the maximum value of rest-mass density at the time (tc = 17413M) when cooling is 
turned on. The asterisk on the curve denotes the value of the central density corresponding to the maximum-mass TOV 
NS (p max, TOV — 2.16 X 10^^ g/cm^). Soon after the maximum density of the TZIO crosses Pmax, TOV, the remnant collapses 
to a BH. Middle panel: Evolution of minimum value of the lapse function (amin) with cooling (solid curve) and without 
cooling (dotted curve). Right panel: Evolution of rest-mass contained within different radii with cooling turned-on. Here 
Mo,r<ro stands for the rest mass contained within a coordinate sphere of radius ro in units of km, centered on the remnant's 
center of mass. These plots demonstrate that the TZIO remnant collapses as a whole. All plots correspond to case Al, and 
M 2.38 Mq 3.52 km 1.17 x 10"^ s. 



the local constraint violations. In addition, the NS (i.e. 
the most compact object in our scenario) remains almost 
unaffected and cold throughout the evolution, i.e., ra- 
diation has no effect on the spacetime structure in the 
vicinity of the NS. In all simulations we present here for 
radii (r) such that r > i^core ^ 20km the local constraint 
violations remain smaller than 0.1% throughout the evo- 
lution after cooling is turned on, justifying our neglecting 
of the radiation self-gravity. 



VI. COOLING OF TZIOS FORMED IN 
HEAD-ON COLLISIONS 

We found in [1] that all our binary pWDNS head- 
on collisions formed hot, quasiequilibrium TZlOs, which 
were more massive than the maximum mass our cold EOS 
can support. However, these remnants did not collapse 
promptly to a black hole. As there is no angular mo- 
mentum involved in a head-on collision, the additional 
support that prevents collapse arises from thermal pres- 
sure alone. Therefore, if one were to cool these objects, 
one would expect that they would eventually collapse on 
a cooling time scale. We check this expectation here so 
that we may implement the same cooling mechanism in 
the inspiral case, where the outcome is not so certain. 

To determine the dominant cooling mechanism, and 
hence the cooling time scale, one needs to know the den- 
sity and temperature of the matter. We estimated the 
temperature of realistic TZlOs to be of order 10^ K. 
Given that typical WD densities are of order 10^ g/cm^, 



it is likely that cooling will be dominated by thermal 
neutrino processes. Realistic neutrino cooling time scales 
are at best of order 1 yr (see discussion in Sec. IVIIip , 
or equivalently ^ 10^ TZIO dynamical time scales. This 
slow cooling rate ensures that the collapse of TZlOs will 
be quasistatic. However, realistic cooling time scales are 
so long that it would be impossible to follow this secular 
phase with hydrodynamic simulations in full GR because 
of the prohibitive computational cost. 

Nevertheless, to confirm that these TZlOs collapse to 
BHs after they have cooled, we can simply scale up the 
cooling law, as long as we keep the cooling time scale 
longer than the hydrodynamical time scale. 

We can estimate the dynamical time scale of a TZIO 

as 

where i^xzio and Mxzio are the radius and mass of the 
remnant, respectively. Here we define the radius of TZlOs 
by the radius of a coordinate sphere that contains 90% 
of the remnant's total mass. For case Al we find txzio ^ 
2300 km ^ 654.5M. The TZIO dynamical time scale in 
case A3 is approximately the same as that of case Al. 
We choose Tc = 6000 km for both cases, so that initially 
Tc ^ 2.6tTZio- This way we reduce the simulation time, 
while maintaining the time-scale inequality. 

Though we choose Tc to be only a few times txzio, 
it is > 100 times larger than the dynamical time scale 
(^core) of the innermost, cold, NS core, which satisfies 
^core/^TZio ~ 15^^^ ^ 60. This property is important 
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because the core of the TZIO collapses first. Finally, 
note that as the entire remnant contracts, its dynamical 
time scale decreases so that the required inequalities are 
better satisfied with increasing time. 

Our expectation is that these TZlOs will collapse fol- 
lowing cooling and that collapse should occur on a cooling 
timescale. Given that (a) the shock-induced thermal en- 
ergy is comparable to the gravitational binding energy, 
and (b) the total remnant mass is close to the maximum 
mass of a cold configuration allowed by our adopted EOS, 
we also expect that a large fraction of the thermal energy 
needs to be removed to induce collapse. Note that this 
expectation applies solely to TZlOs formed in head-on 
collisions. When TZlOs form after the merger of WDNSs 
in an initially circular orbit, the remnant rotates. Hence, 
we cannot tell a priori that collapse will take place be- 
cause of the additional centrifugal support. 

Table |II] summarizes the physical parameters of the 
head-on collision cases we study, and Table IIIII presents 
the AMR grid structure used in each case. All our simu- 
lations with cooling turned on show that the TZlOs rem- 
nants formed in the head-on collisions begin to contract, 
and within a few cooling time scales collapse to a black 
hole. In contrast, we find that when cooling is turned off 
the remnant does not collapse and remains in quasiequi- 
librium. All these results can be clearly seen in Figs. [1] 
and [21 which correspond to case Al. Case A3 is qualita- 
tively similar. 

In Fig. [1] rest-mass density contours are plotted in the 
orbital plane at selected times. Figure [2] shows the evolu- 
tion of the maximum rest-mass density with and without 
cooling, the minimum value of the lapse function, and the 
rest-mass contained within different radii from the center 
of mass. If the cooling mechanism remains off, both the 
maximum value of the density and the minimum value of 
the lapse remain constant with time (see left and middle 
panels of Fig. [2]). We find that the same holds true for 
the rest-mass contained within 40 km and 220 km. 

By contrast, if cooling is turned on, the maximum 
density (minimum lapse) increases (decreases) with time. 
Moreover, the rest mass contained within 40 km and 220 
km increases with time, indicating that the outer layers 
are also contracting as the TZIO cools (see right panel 
in Fig. [2]). Initially, the maximum density (minimum 
lapse) increases (decreases) almost linearly with time, un- 
til Po,max crosses the value of 2.16 x 10^^ g/cm , which 
corresponds to that of a maximum mass NS configura- 
tion built with our cold EOS. Soon after this point, the 
remnant essentially free-falls, the density blows up, the 
lapse function plummets and a BH is eventually formed. 
The BH in case Al can be seen in Fig. [3l where we plot 
rest-mass density contours and K = P/Pcoid contours in 
the orbital plane in the innermost 12 km of the remnant, 
which contain about 85% of the total mass at the time 
of BH formation. The K contours show that the matter 
around the BH is cold, i.e., ~ 1, as expected. 

Cases Al and A2 collapse to a black hole after about 
5 and 3 coohng time scales, respectively, which is ex- 



pected as the collapse proceeds without additional shock 
heating. The mass of the black hole when an apparent 
horizon forms for the first time is Mbh ^1.4 Mq and 
the coordinate radius of the BH (in our adopted gauge) 
is Rbu ^ 1-08 km ^ 0.53Mbh. Thus we have demon- 
strated that our cooling mechanism yields results which 
are consistent with our theoretical expectations. 



VII. BINARY WDNS INSPIRAL 

To predict the final outcome of a binary WDNS in an 
initially circular orbit, we performed a simulation of a 
corotating binary pWDNS starting at the Roche limit 
separation. Throughout, we label this case by the letter 
A. Table [TTl outlines the physical parameters of case A, 
and Table Hill outlines the adopted AMR grid structure. 

For the simulations performed here, we were able to 
demonstrate 2nd-order convergence for the first quarter 
of an orbit monitoring the conservation of angular mo- 
mentum, and the constraint violations. The convergence 
study showed that angular momentum decays linearly 
with time, but the linear decay rate decreases with in- 
creasing resolution to second order. Moreover, this de- 
cay rate remains roughly constant up until merger. Fur- 
thermore, a resolution study using pWDNSs systems has 
been carried out in [1] , where we showed that the results 
were qualitatively insensitive to resolution implying that 
the resolutions used were sufficiently high. The resolu- 
tion used in our inspiral pWDNS calculations is twice 
that used in [1] indicating that our simulations are well 
within the convergent regime. 

A. Initial configuration 

We prepared valid general relativistic initial data as 
described in Section IIIII The ADM masses of the com- 
pact objects in isolation we consider are 0.98 Mq and 
1.4 Mq for the pWD and NS, respectively. After solving 
the CTS equations, we map po, ^ and a, /3\ and v\ from 
the grids used in the elliptic solver code onto the grids 
used in the evolution code via second-order polynomial 
interpolation. For accuracy, we make sure that the reso- 
lution of the initial data grids is always higher than the 
resolution of the evolution grids. 



B. Dynamics of the WDNS merger 

In [s^l we analyzed the stability of corotating binary 
WDNSs at the Roche limit, accounting for GR effects on 
the mass-radius relationship of the WD. We concluded 
that if the mass ratio q = Mwd/^ns is larger than a 
critical mass ratio g'crit ^ 0-5, then mass transfer from 
the WD to the NS will be unstable, and the WD will 
be tidally disrupted. The binary pWDNS system simu- 
lated in this work has a mass ratio = 0.7, so we expect 
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FIG. 3. Left: Snapshot of rest-mass density profile at the end of the Al simulation shown in Fig. (TJ Contours represent the 
rest-mass density in the orbital plane, plotted according to po = po,maxlO~°'^-^~°"°^^ (j = 0, 1, . . . , 7). Here y9o,max = 4.645pnuc, 
where pnuc = 2 x lO^'^ g/cm^. Right: Snapshot of K = P/Pcoid profile at the end of the Al simulation. The contours represent 
K in the orbital plane, plotted according to K = ]^Q-o.ii25j+o.9 _ 1, . . . , 7). The plot demonstrates that the matter near 
the BH is cold {K ^ 1), as expected, and K increases with the distance from the core. The color coding is the same as that 
used in Fig. [T] with white indicating K ^ 1. The black disk in the center denotes the BH apparent horizon. The plots focus 
in the innermost 12 km from the TZIO center of mass, where the object is approximately spherical, and for this reason we do 



not show XZ and YZ meridional slices, 
isolated stars. 



Here M = 2.38 M© = 3.52 km = 1.17 x 10"^ s is the sum of the ADM masses of the 



that the system should experience tidal disruption and 
merge on an orbital time scale soon after mass transfer 
has started. 

In our simulations, the pWDNS binary completes al- 
most 2.5 orbits before the pWD is completely disrupted. 
Figure |4] plots rest-mass density contours in the orbital 
plane at selected times for case A. The top row, middle 
panel shows the binary shortly after completing the first 
orbit. At this time, an accretion stream from the pWD 
to the NS develops, followed by the formation of an ac- 
cretion disk around the NS. Matter from the accretion 
stream smashes into the accretion disk, shock heating 
the gas at that location. This process continues until the 
pWD is completely disrupted. After pWD tidal disrup- 
tion, a long tail forms that moves outwards and around 
the NS. The pWD matter that orbits the NS at closer 
separations collides with the tail inducing further strong 
shocks. 

The bottom row, left panel of Fig. |4] shows the system 
after 3 orbits have been completed. At this point, the 
pWD is completely disrupted and a large, rotating, mas- 
sive mantle and disk around the NS has begun to form. 
A tidal tail is also visible. This snapshot is followed by a 
long epoch in which the rotating mantle settles onto an 
extended disk around the central object, composed of a 
slowly spinning, cold NS core surrounded by a hot atmo- 



sphere and disk composed of pWD debris. We find that 
even at this late stage, the NS core maintains its original 
spin, and the hot mantle surrounding it spins and settles 
into quasiequilibrium. Non-axisymmetric clumps of mat- 
ter inspiraling near the cold NS core launch spiral density 
waves into the disk. The remnant of the pWDNS merger 
may be best characterized as a spinning TZIO with an 
extended (i^disk ~ 1000 km), massive disk. 

We define the radius of the TZIO (i^xzio) as the dis- 
tance between the center of mass and the "north" pole 
of the remnant. The z-radius of the remnant for a cut- 
off density ~ 10~^po,max, where po,max is the maximum 
density of the remnant, is i^xzio ^ 300 km. We esti- 
mate the rest mass of the TZIO (Mxzio) as the rest mass 
contained within a sphere of coordinate radius equal to 
^TZio, and the disk mass (Mdisk) by subtracting Mxzio 
from the total rest mass. We find Mxzio = 1.82 Mq and 
Mdisk = 0.7 Mq. The disk is massive and > 50% of the 
original WD rest mass is eventually stored in the disk. 

To first order, the TZIO is spherical. This is evidenced 
by the XY, XZ, and YZ, density contour plots of Fig. [5l 
which focus on the innermost regions of the remnant at 
the end of the simulation. 

The cutoff density in all the density contour plots we 
show here is 10~^"^ km~^ ^ 10~^'^po,max, which is ap- 
proximately 4 decades above atmosphere density. The 
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FIG. 4. Snapshots of rest-mass density profiles at selected times for case A. The contours represent the rest-mass density in the 
orbital plane, plotted according to po — y9o,maxlO~°'^^'^~°"^^ (j = 0, 1, . . . , 9). The insets focus on the NS and demonstrate it is 
nearly unaffected by the merger. The inset density contours are plotted according to po = Po,maxl0~'^'^^^'^~°"^^^ (j = 0, 1, . . . , 7) 
In both cases Po,max — 4.645/?nuc, where />nuc = 2 x 10^^ g/cw? . The color coding is the same as that used in Fig.[Tl with white 
indicatine: near vacuum. Here M 2.38 M© 3.52 km 1.17 x 10"^ s is the sum of the ADM masses of the isolated stars. 



equatorial and polar coordinate radii of these contours 
is Ve ^ 350 km and Vp ~ 150 km, respectively. There- 
fore, the ratio of these radii is Vp/ve ~ 3/7. Given that 
the core of the remnant is approximately spherical, the 
smallness of the ratio r^/rg indicates that the disk stores 
a large amount of angular momentum. 

In the bottom row of Fig. [5l we plot contours of 
K = P/Pcoid- These entropy contours show that the 
neutron star core is cold K ^ 1 and that K increases 
with the distance from the core in the orbital plane and 
with increasing z in meridional planes. This is reminis- 
cent of the K pattern observed in our binary pWDNS 
head-on collision studies in [1], where K ^ 1 in the core, 
but increases with distance from the center. Note the 
spiral density wave pattern visible in the bottom row, 
left panel of Fig. [5l 

Unlike the head-on collision case, in which the out- 
ermost layers of the NS are shock heated and stripped 
away when the NS smashes into the denser parts of the 
pWD (Fig. 4 insets, [T]), after a pWDNS binary inspiral 
and merger, the NS retains its outer layers, and its struc- 



ture remains nearly unaffected throughout the simulation 
(Fig. H] insets). Moreover, in the head-on collision, about 
18% of the total initial rest mass is ejected to infinity, 
but in the inspiral case, no ejection of matter to infinity 
is observed. 

As in the case of head-on collisions, we find that the 
typical temperature in our inspiraling binary remnant is 
of order 10^^ K. In Fig.[6]we show temperature profiles of 
the remnant. To estimate the temperature T, we assume 
that the temperature dependence of eth can be modeled 
as 



3kT 
2mn 



J 7 

po 



(68) 



where rrin is the mass of a nucleon, k is Boltzmann's con- 
stant and a is the radiation constant. The first term rep- 
resents the approximate thermal energy of the nucleons, 
and the second term accounts for the thermal energy due 
to relativistic particles. The factor / reflects the number 
of species of relativistic particles that contribute to the 
thermal energy. When T <C 2me/k ~ lO^^K, where rrie 
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FIG. 5. First row: Snapshots of rest-mass density profiles at selected times for case A. The contours represent the rest- 
mass density in the orbital plane and the XZ and YZ meridional planes, plotted according to po = />o,maxlO~°'^^-^~°"^^ (j = 
0, 1, . . . , 7), where po,max = 4.645pnuc, and pnuc = 2 x 10^^ g/cm^. Second row: Snapshots oi K = P/Pcoid profiles at selected 
times for case A. The contours represent K in the orbital plane and the XZ and YZ meridional planes, plotted according to 
K = io-o-ii25i+o.9 ^ 0, 1, . . . , 7). The plots show that the remnant NS core is approximately spherical and cold {K ^ 1). 
Far from the core the remnant is hot. K increases as we move away from the core and the orbital plane. All plots focus in the 
innermost 200 km from the TZIO center of mass. The color code used is the same as that defined in Fig. [1] with white color in 
the second row indicating K ^ 1. Here M = 2.38 M© = 3.52 km = 1.17 x 10"^ s. 



is the electron mass, thermal radiation is dominated by 
photons and / = 1. When T ^ 2me//c, electrons and 
positrons become relativistic and also contribute to radi- 
ation, and / = 1 + 2 X (7/8) = 11/4. At sufficiently 
high temperatures and densities (T > 10^^ K,po ^ 
lO^^g cm~^), neutrinos are generated copiously and be- 
come trapped. So, taking into account three flavors of 



neutrinos and antineutrinos, / = 11/4 + 6 x (7/8) = 8. 
In our temperature estimate, we find / self-consistently 
in the following sense: we first calculate the temperature 
assuming / = 0. If the calculated temperature and den- 
sity are inconsistent with our choice of / (which we test 
based on the above inequalities), then we choose a dif- 
ferent /, until we find the value of / which is consistent. 
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TZlOs formed in head-on cohisions (see Section |VI]) , and 
compare the result to an uncooled case. 
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FIG. 6. Temperature (T) profile for case A. The temperature 
is presented in units of 10^^ K, where K indicates degrees 
Kelvin (not to be confused with the EOS entropy parameter 
K) . The solid curve corresponds to the values of T at the end 
of the simulations and along the x-axis for ^ = ^c, x > Xc, 
where Xc (ye) is the x (y) position of the center of mass of 
the remnant. The dotted curve corresponds to the values of 
T along the ^/-axis for y > yc^ x — Xc. It is clear that typical 
temperatures are of order 10^^ K, while the TZIO core is 
practically at OK. For a realistic massive WDNS merger we 
expect T ^ 10^ K (see discussion following Eq. ()71|) ). 



We find that / = 11/4 is consistent with the tempera- 
tures and densities in our pWDNS merger. However, we 
expect that in realistic mergers / = 1, as the expected 
temperatures are of order l{fK (see discussion following 
Eq. m)- 

To solve Eq. ()68p for T we need to know eth- We cal- 
culate eth via 



eth 



{K - l)Peold 

(Fth -l)po'- 



(69) 



where Eqs. (|37|) . (|38|) and the definition of K were used 
to obtain this last equation. 

As in our pWDNS head-on collision studies, we find 
that the pWDNS binary remnant does not collapse 
promptly to a BH. In contrast to the head-on collision 
cases, in which only thermal pressure supported the rem- 
nant from prompt collapse, collapse may be delayed in 
the pWDNS binary inspiral case by both thermal pres- 
sure and centrifugal support. One way to assess the im- 
portance of thermal support in the binary inspiral rem- 
nant is to apply the same cooling recipe as in the case of 



Cooling of spinning TZlOs formed in WDNS 
mergers 



To determine whether the spinning TZIO from the 
WDNS inspiral and merger will collapse to a BH fol- 
lowing coohng, we apply our cooling technique setting 
the cooling time scale to the same value used earlier 
(tc = 6000 km) in Sec. IVII We turn on cooling after 
about 4.5 orbital time scales and follow the subsequent 
evolution for about 7 cooling time scales. 

Figure [71 plots density contours in the orbital plane 
at selected times, and Fig. [8] shows the evolution of the 
maximum rest-mass density and the rest mass contained 
within spheres with different coordinate radii. These fig- 
ures demonstrate that both the maximum density, and 
the rest mass within 40 km and 220 km increase as func- 
tions of time. As the hot outermost layers become cooler 
they contract and accumulate onto the innermost colder 
parts. For this reason the innermost density contours of 
like density increase in size (see Fig. [7]). The remnant is 
contracting with time. The contraction is more rapid in 
the beginning and begins to plateau after about 6 cool- 
ing time scales. Therefore, the spinning TZIO does not 
collapse to a BH when thermal support is removed. 

Figure [9] shows density and K contours for the inner- 
most regions of the remnant in the XY, XZ, and YZ 
planes, 6 cooling time scales after cooling was turned on. 
The shape of the NS core remains spherical throughout 
the evolution. The XZ and YZ contours demonstrate 
that the disk and mantle have become thinner, as ex- 
pected when cooling takes place. Due to this effect, this 
final configuration is massive accretion disk onto a NS, 
rather than a disk around a TZIO. 

The bottom row of Fig. [9] shows contours oi K = 
P/Pco\d' Notice that the neutron star core remains cold 

^ 1 at the end of the simulation and that elsewhere 
K has decreased considerably compared to the run with- 
out coohng. Here iiTmax ^ 1-25, while in the run without 
cooling i^max ^ 10. In the innermost region, cold pres- 
sure dominates, with K ^ 1.05. Given that the rest mass 
within 220 km of the remnant center of mass is greater 
than 2.05 M©, which exceeds the maximum supportable 
mass by our cold EOS, we conclude that the spinning 
TZIO is centrifugally supported from collapse to a BH. 

Based on these results and the scalability of our simu- 
lations to the reahstic scenario, we are led to the tentative 
prediction that realistic WDNS mergers with total rest 
mass < 2.5M0, the rest mass in our simulations, will not 
collapse to a BH following cooling. This conclusion as- 
sumes that angular momentum redistribution takes place 
on a longer time scale than cooling. 

Given the absence of outfiows in our simulations (with 
the caveat that we do not model nuclear reactions), the 
final total rest mass (?^ 2.5 Mq) is larger than the maxi- 
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FIG. 7. Snapshots of rest- mass density profiles at selected times for case A with cooling turned on. The contours represent the 
rest-mass density in the orbital plane, plotted according to po — po,maxlO~°'^-^~°'^^ (j = 0, 1, . . . , 8), where po,max = 4.645pnuc, 
and yOnuc = 2 X 10^^ g/cm^. The last two snapshots show that the contraction practically stops after about 6 cooling time 
scales. The color code used is the same as that defined in Fig. [T] and M — 2.38 Mq — 3.52 km = 1.17 x 10~^ s. 



mum rest mass supportable by our cold EOS, even allow- 
ing for maximal uniform rotation. Therefore, we expect 
that after viscosity and/or magnetic fields redistribute 
angular momentum, the remnant will collapse to a black 
hole. This conclusion will be true in the case of realistic 
WDNS mergers, unless the true nuclear EOS supports a 
uniformly rotating star with a rest mass exceeding the 
remnant mass. Many viable EOSs do not support a uni- 
formly rotating cold configuration with rest mass as large 
as 2.5M0 the remnant rest mass in our simulations. 



VIII. DISCUSSION 



To identify the relevant nuclear reaction networks and 
the dominant cooling mechanisms in realistic inspiraling 
WDNS binaries, we need to estimate the temperatures of 
realistic TZlOs. Moreover, to determine the time scale 
on which angular momentum redistribution occurs we 
have to consider viscosity and/or magnetic fields. In this 
section we discuss these issues. 
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A. Temperature 

The characteristic temperature of reahstic TZlOs is 
expected to be of order 10^ K. This is because the en- 
ergy available for shock heating is of order the gravita- 
tional interaction energy when the two stars first touch, 
Mns^wd/^wd- Our simulations demonstrate that the 
NS is largely unaffected by shock heating and remains 
cold. Hence, most of the thermal energy is stored in the 
WD debris. The total thermal energy, £^th, is then 



Eth 



MwD 



MnsMwd 



(70) 



From this last equation we can estimate the characteristic 
temperature as 



qk 



(71) 



where CwD = ^wd/^wd is the WD compaction. All 
things being equal (i.e., no mass loss, same mass ratio, 
etc.), characteristic TZIO temperatures should be pro- 
portional to the WD compaction. The compaction of a 
realistic 1.0 Mq WD that obeys the Chandrasekhar EOS 
is CwD ^ 10-"^ [34, 68]. If the NS mass is 1.4 M© then 
q 0.7, and Eq. predicts T 1.55 x 10^ K. 

Note that applying Eq. ([TT]) to case A, where CpWD — 
10~^, yields a temperature T 1.55 x 10^^ K, i.e., in 
good agreement with our simulations [89i] . 



B. Nuclear fusion 

Are realistic WDNS binary remnant densities and tem- 
peratures high enough for nuclear reactions to take place? 
The shock-heated matter is composed of hot (diluted) 
WD debris, so its density is of order typical WD densi- 
ties, i.e., 10^ g/cm^. A 1.0 Mq WD is sufficiently massive 
that its main constituent elements are carbon and oxy- 
gen. While the temperatures and densities we expect for 
realistic mergers are probably not high enough for oxy- 
gen burning to become important, they are sufficiently 
high for carbon fusion to become dynamically relevant. 

Non-explosive nuclear reactions in the context of 
WDNS mergers were recently considered in [90]. A ID 
steady state model of accretion onto a NS was introduced, 
allowing for disk wind outflows that do not exert any 
torque on the disk. It was found that heating from nu- 
clear burning is so important that a disk wind eventually 
unbinds 50% — 80% of the original WD mass. It was 
suggested that these eject a may include small quantities 
of radioactive ^^Ni. In such scenarios, detectable EM 
signals will likely follow a WDNS merger. Although this 
ID steady-state model includes much of the important 
physics (albeit in parametrized form), it is simplified and 
does not apply to the large mass-ratio mergers simulated 
here. However, as in the ID model we do expect that 
nuclear burning will also be non-explosive in a realistic 
WDNS merger, as we now explain. 



In a head-on collision of a WDNS binary with compan- 
ions of comparable mass that collide at free-fall velocity, 
the kinetic energy of motion is converted by shocks into 
thermal energy in the WD remnant. This shock heating 
at merger guarantees that a degenerate WD initially in 
hydrostatic equilibrium will acquire shock-induced ther- 
mal pressure comparable in magnitude to its original 
equilibrium degeneracy pressure, thereby lifting the de- 
generacy, i.e.. 
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(72) 



The net effect should be to reduce the likelihood of ex- 
plosive carbon burning, since a carbon flash requires a 
degenerate environment. The reason for this is that if 
gas pressure becomes a significant component of the to- 
tal pressure, then the pressure will be sensitive to the 
temperature. Therefore, if carbon fusion takes place, 
the released heat will increase the gas temperature which 
will, in turn, increase the pressure. As a result the gas 
will expand, decreasing its density and temperature, and 
eventually carbon fusion will be turned off. Such a pro- 
cess is self-regulated, a well-known fact. 

Shock heating plays a similar role in the merger of an 
inspiraling binary, only it is not as strong and the frac- 
tion of the thermal pressure generated will be smaller, 
due to the role of angular momentum in lessening the 
impact and contributing to the support of the remnant. 
Using our estimated temperature T ^ 10^ K and charac- 
teristic density 10^ g/cm^ for realistic TZlOs, the ratio 
of thermal gas pressure to the electron degenerate pres- 
sure is 4. This implies that the WD debris would be 
non-degenerate. Under these conditions a carbon flash is 
likely suppressed, but further simulations would be useful 
to confirm this. 



C. Neutrino cooling 

For the estimated characteristic temperature T = 
10^ K and density po = 10^ g/cm^, the dominant cool- 
ing mechanism likely will involve neutrino emission. At 
these densities and temperatures, thermal neutrino pro- 
cesses (pair neutrinos, photoneutrinos, plasmon decay, 
and bremsstrahlung [68]) are important, with pair anni- 
hilation (e~ +e+ — ^ being slightlymore important 
than the other processes (see Fig. 1 in [91| and Fig. 3(a) 
in [921). The pair neutrino cooling rate can be estimated 
as [93] 



^pair ^ 4 45 ^ lO^I^ [erg/g/s], (73) 

P6 

where Tg = T/10^ K, pe = po/lO^g cm-^ and the high- 
temperature and non-degenerate limit has been assumed. 

The specific thermal energy is approximately given by 
eth = SkT/rrin' Based on this, the cooling time scale can 
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FIG. 8. Left: Evolution of maximum value of rest-mass density with cooling (solid curve) and without cooling (dotted curve) 
for case A. Here y9o,max is the maximum value of the rest-mass density, />o,max (te) = 5.88 X 10-^ km-2 = 7.92 x 10^^ Qi/cm^ is 
the maximum value of rest-mass density at the time (tc = 36705M) when cooling is turned on. Right: Evolution of rest-mass 
within different radii with cooling turned-on. Here Mo,r<ro stands for the rest mass contained within a coordinate sphere of 
radius ro in units of km, centered on the remnant's center of mass. These plots demonstrate that the TZIO remnant contracts 
when cooling is turned on. However, the values of po,max and Mo,r<ro begin to plateau after 6 cooling time scales, indicating 
no further contraction proceeds after this time. Here M = 2.38 Mq = 3.52 km = 1.17 x 10~^ s is the sum of the ADM masses 
of the isolated stars. 



be estimated as 

'^cooli] 



pair 
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(74) 



Notice that for Tg = 0.1, Tcooiing ^ 10^ yr. Thus the 
object cools fast when it is very hot, but when the tem- 
perature drops to lO^K it takes hundreds of millions of 
years for cooling to take place. 

Based on these considerations and Eq. ([73]), the net 
conclusion is that the neutrino cooling time scale is highly 
temperature sensitive, and our pWDNS inspiral simula- 
tion may only provide a crude estimate of temperature. 
Therefore, simulations with more physics are necessary to 
precisely calculate realistic TZIO temperatures, so that 



the relevant cooling time scales may be better estimated. 

Adopting the cooling rate ([73]) we can estimate 
whether neutrinos from WDNS mergers are detectable. 
The number of detectable neutrinos (Nd) are approxi- 
mately given by 



(75) 



where Lj^ is the total neutrino luminosity, the neu- 
trino detection cross section, AT the time interval over 
which neutrinos are emitted, D the distance to the bi- 
nary, and e^y the average neutrino energy. Given that 
^ 10~^^cm^, and the expected neutrino energy from 
pair annihilation is ~ 0.5MeV, we estimate 
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(76) 



where we have assumed that the entire TZIO mantle 
emits neutrinos at the same rate for a year. 

Given this result, we conclude that neutrinos emitted 
in WDNS mergers are unlikely to be detectable. How- 
ever, simulations with detailed microphysics [94] would 
be useful to confirm this. 



D. Angular momentum redistribution 



Our inspiraling WDNS merger simulation with cool- 
ing turned on shows that the remnant does not collapse 
to a BH following cooling, because it is centrifugally sup- 
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FIG. 9. First row: Snapshots of rest-mass density profiles at selected times for case A with cooling. The contours represent the 
rest-mass density in the orbital plane and the XZ and YZ meridional planes, plotted according to po = pcmaxlO"^"^^"'"^"^^ (j — 
0, 1, . . . , 7), where y9o,max = 4.645pnuc, and pnuc = 2 x 10^^ g/cm^. Second row: Snapshots oi K — P/Pcoid profiles at selected 
times for case A with cooling. The contours represent K in the orbital plane and the XZ and YZ meridional planes, plotted 
according to = ^^q-o.omj+o.i (j ^ q, 1, . . . , 7). The plots show that the remnant NS core is approximately spherical and cold 
[K ^ 1). Far from the core the remnant is hotter. K increases as we move away from the core and the orbital plane. In contrast 
to the case without cooling (see Fig. O, the maximum value for K here is i^max ~ 1.25. For easy comparison with the case 
without cooling (Fig.[5|), all plots focus in the innermost 200 km from the remnant center of mass. The color code used is the same 
as that defined in Fig.[Tl with white color in the second row indicating K ^ 1. Here M = 2.38 Mq = 3.52 km = 1.17 x 10~^ s. 



ported. Given that the mass of the remnant is larger than 
the maximum mass supportable by our cold EOS, it is 
likely that delayed collapse will take place after angular 
momentum is redistributed. 

Angular momentum redistribution will occur on the 
viscous or Alfven time scale. Assuming an a-disk, the vis- 



cous time scale (neglecting the disk self-gravity) is given 

by 



i?3 



TZIO 



(77) 
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where H is the disk scale height, R the characteristic disk 
radius, and a the turbulent viscosity parameter. Using 



the values for H/ and Mxzio found in our simulations 
we estimate that in realistic WDNS mergers the viscous 
time scale is 



R y^V MTzio 



1/2 



(78) 



where R ~ 2i?wD, i-e., near the Roche limit for a l.OM© 
WD with a lAM^ NS. 

The Alfven time scale {Ia = R/va^ where va is the 
Alfven speed) is given by 



i?3 



Mtzio ' 



where the /3 parameter 



8nP 



(79) 



(80) 



was introduced to obtain the last expression. If we use 
the same values for i?, H/R and Mxzio as in Eq. (|78]) . 
the Alfven time scale becomes 



R 



Mtzio ^ -^/^ 

104 km/ Vl-8AfQ 



(81) 



We emphasize that the dimensionless parameters a and 
P above are unknown and may both be <C 1, in which 
case the angular momentum redistribution time scale 
may be as long as the cooling time scale. For exam- 
ple, observations of magnetic WDs indicate that surface 
magnetic field strengths are B ~ 10^ — 10^ G [95] or 
P ~ 10~^^ — 10~^ <C 1, where we calculated the ther- 
mal pressure as P = pokT/rrin with po = 10^ g/cm^, 
T = 10^ K. If /3 < 10-^^ then the Alfven time scale is 
longer than 1 yr, i.e., the cooling time scale for Tg = 1. 
However, field amplification via winding and instabilities 
(e.g. magnetorotational instability) is always possible. 
Hence, we must await detailed calculations for reliable 
estimates of the angular momentum redistribution time 
scale. 



IX. SUMMARY AND CONCLUSIONS 

This work is a follow-up to our study of binary WDNS 
head-on collisions [1], focusing on the dynamics of an ini- 
tially circular, quasiequilibrium WDNS binary through 
inspiral and merger. In particular, we begin with a cir- 
cular binary in which the WD has just filled its Roche 
lobe (the Roche limit) and with systems whose total mass 
exceeds the maximum mass that a cold EOS can support. 
The goal is to determine whether a WDNS merger leads 
to either prompt collapse to a BH or a spinning quasiequi- 



I 

librium configuration consisting of a cold NS surrounded 
by a hot gaseous mantle of WD debris, or something else. 

Due to the vast range of dynamical time and length 
scales, hydrodynamic simulations in full GR of realistic 
WDNS mergers (head-on or otherwise) are computation- 
ally prohibitive. For this reason, we tackle the problem 
using the same approach as in our investigation of bi- 
nary WDNS head-on collisions. In particular, we adopt 
the pseudo-white dwarf (pWD) approximation with the 
10:1 EOS constructed in [l|. This EOS captures the main 
physical features of NSs, but scales down the size of WDs 
so that the ratio of the isotropic radius of a TOV 0.98 Mq 
pWD to that of a TOV 1.5 Mq NS is 10:1 (hence the 
name of the EOS), rather than the more realistic ra- 
tio 500:1. These pWDs enable us to reduce the range 
of length and time scales involved while maintaining all 
length and time-scale inequalities, rendering the compu- 
tations tractable and the results scalable. 

If the pWDNS merger does not result in prompt col- 
lapse to a black hole, it is unlikely that the corresponding 
WDNS merger will collapse promptly. 

The reason for this expectation is that the pWD ap- 
proximation is based on scaling. In particular, both the 
colHsion velocity and the pre-shocked WD sound speed 
scale as ~ (M/i^wo)^^^- This implies that the Mach 
number is invariant under scaling of Rwb and so is the 
degree of shock heating. So the thermal energy, as well 
as the rotational kinetic energy (T) and the gravitational 
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potential energy (W) all scale as ~ M'^/Ry^/i:), when the 
binary merges. Thus T/|W^| is also invariant under scal- 
ing of RwD • These considerations simply mean that with 
respect to gravity the relative importance of thermal and 
rotational support in a WDNS merger remnant is approx- 
imately invariant, when the masses of the binary compo- 
nents are fixed and the only quantity that changes is the 
WD radius. As a consequence, the results obtained when 
adopting pWDNS systems can be scaled up to realistic 
WDNS systems. 

To predict whether a TZIO, which does not collapse 
to a BH promptly, will collapse following cooling, we in- 
troduced an artificial cooling mechanism (see Sec.|V]). If 
following cooling the remnant collapses, we expect that 
delayed collapse in the corresponding WDNS case likely 
will take place on a cooling time scale. 

To test our cooling prescription, we applied it to the 
TZlOs formed in the WDNS head-on collision simulations 
we performed in py. We demonstrated that these rem- 
nants collapse to a black hole when the excess thermal 
energy is radiated away, as expected. 

Finally, we simulated the merger of an initially 
quasiequilibrium, corotational pWDNS system in circu- 
lar orbit at the Roche limit, comprised of a IAMq NS 
and a O.OSM© pWD. We find that the remnant of the 
pWDNS inspiral is a spinning TZIO which is surrounded 
by an extended, hot disk. The coordinate radius of the 
TZIO remnant and disk is approximately 300 km and 
1000 km, respectively. We estimated the disk mass to be 
> 50% of the initial original WD rest mass. In contrast 
to our binary WDNS head-on collision investigations, no 
outfiows were observed in the circular case. The final 
total ADM mass (~ 2.4 Mq) is greater than the maxi- 
mum mass supportable by a cold, degenerate star with 
our adopted NS EOS. However, the remnant does not 
collapse promptly to a black hole. This is because the 
remnant is both thermally and centrifugally supported. 
To determine whether centrifugal support by itself sup- 
ports the remnant from collapse, we enabled our radia- 
tive cooling mechanism and found that the object does 
not collapse to a black hole following cooling. Therefore, 
the extra support provided by rotation is sufficient for 
holding the collapse. 

Although the TZIO does not collapse following cool- 
ing, ultimate collapse to a BH is almost certain, since 
the final total mass is larger than the maximum possi- 
ble mass supportable by our cold EOS (and many nu- 
clear EOSs), even allowing for maximal uniform rota- 



tion. Therefore, delayed collapse likely will take place 
after viscosity or magnetic fields redistribute the angu- 
lar momentum and/or following cooling. This conclusion 
will be true in the case of realistic WDNS mergers, unless 
the true nuclear EOS supports a uniformly rotating star 
with a rest mass exceeding the remnant mass. Many vi- 
able EOSs do not support rest masses as large as 2.5Mq 
[37], the remnant rest mass in our simulations. 

Our results hold true provided that nuclear burning 
remains unimportant in the post-merger event. We esti- 
mated that typical realistic TZIO temperatures will be of 
order 10^ K. For typical WD densities of order 10^ g/cm^ 
carbon is ignited and can become an important source of 
heating. Though nuclear burning likely will play some 
role in the post-merger evolution of a massive WDNS 
system, we do not expect a carbon fiash to occur. The 
reason for this is that shocks at merger lift the degeneracy 
of the WD matter. Given that a carbon fiash requires a 
cold degenerate environment, the net effect should be to 
reduce the likelihood of explosive carbon burning. Nev- 
ertheless, further simulations would be useful. 

The neutrino cooling time scale is highly temperature 
sensitive, and our pWDNS inspiral simulation may only 
provide a crude estimate of temperature. Therefore, sim- 
ulations with more physics are necessary to precisely cal- 
culate realistic TZIO temperatures, so that the relevant 
cooling time scales may be better determined. Finally, 
while our simulations indicate that prompt collapse to a 
black-hole is not possible for WDNS systems with total 
rest mass < 2.5M0, it is likely that systems with greater 
mass can collapse promptly. Therefore, more simulations 
in full GR are necessary before a definitive solution to the 
problem can be given. We plan to address these issues in 
a future work. 
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